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. Abstract 

P_f The planet Kepler-lQb is known to follow a circumbinary orbit around a 

double system of two main- sequence stars. We construct stability diagrams in 
the "pericentric distance — eccentricity" plane, which show that Kepler-16b 
is in a hazardous vicinity to the chaos domain — just between the instability 
"teeth" in the space of orbital parameters. Kepler-lQb survives, because it 
is close to the half- integer 11/2 orbital resonance with the central binary. 
The neighbouring resonance cells are vacant, because they are "purged" by 
Kepler-lQb, due to overlap of first-order resonances with the planet. 
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A circumbinary planet in the Kepler-16 system was discovered in 2011 
basing on the data from the Kepler spacecraft pQ. Its orbital parameters were 
determined in [TJ [2] . An application of an empirical numerical-experimental 
criterion [3] of the stability of planetary circumbinary orbits, as well as long- 
term numerical integrations of the Kepler-lQb orbit, identifies this system as 
stable, though not far (in about 10-20% in the semimajor axis of the planet) 
from the inner instability zone [HE]- However, a fine structure of the stability 
border in the space of orbital parameters has not yet been studied. 

We adopt the values of the orbital parameters of the double star Kepler- 
16 A-B as given in [TJ, and compute the outer planetary orbits in the planar 
elliptic restricted three-body problem. At the initial moment of time, the 
relative location of the three bodies is that determined in pQ for a particular 
reference epoch. We vary the planet orbital eccentricity e and pericentric 
distance q = a(l — e) (where a is the semimajor axis of the orbit) and 
integrate test planetary orbits on a fine grid of the starting values of q and 
e. 

To explore the stability of the planetary motion, we use two stability 
criteria. The first criterion is the value of the maximum Lyapunov exponent. 
The second criterion is the "escape-collision" one: the orbit is stable if the 
distance between the planet and one of the stars does not become less than 
10~ 3 AU or does not exceed 10 3 AU. 

The Lyapunov exponents of a trajectory measure the rate of exponential 
divergence of nearby trajectories in the phase space of motion [TT] . Let 
d{to) <C 1 be the initial displacement of a "shadow" trajectory from the 
"guiding" one, and d(t) be the displacement at time t. Then the Lyapunov 
exponent is defined by the formula 

L= i im ^_i n 4^L (1) 

Depending on the direction of the initial displacement in the phase space, the 
quantity L of a trajectory of a Hamiltonian system can attain 2N generally 
different values, where N is the number of degrees of freedom. These values 
form symmetric pairs: for each Lj > there exists L i+ ^ = —Li < 0; % — 1, 
. . . , N [11] . Therefore, in practice it is sufficient to compute solely N expo- 
nents Li > 0. The set of all 2N exponents is called the spectrum of Lyapunov 
exponents, or, the Lyapunov spectrum. A non-zero value of the maximum 
(in the spectrum) Lyapunov exponent indicates chaotic character of the mo- 
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tion, whereas zero value indicates that the motion is regular (quasiperiodic 
or periodic) [T2T ITT] . On an everywhere dense set of starting data for shadow 
trajectories, the Lyapunov exponent attains the single (maximum) value, - 
the maximum Lyapunov exponent, which we denote as L in what follows. 
The inverse of this quantity, Tl = is the so-called Lyapunov time. It 
represents the characteristic time of predictable dynamics. 

For computing the Lyapunov spectra we use the algorithms and software 
developed in [TBI [HJ [15] on the basis of the HQRB method by von Bremen et 
al. [16]. This method is based on the QR decomposition of the tangent map 
matrix using the Householder transformation. The trajectories are integrated 
using the integrator by Hairer et al. [T7]. It is an explicit 8th order Runge- 
Kutta method due to Dormand and Prince, with the step size control. 

For separating the regular and chaotic orbits of a dynamical system, a 
statistical method was proposed in [HI [191 120]. It consists of four steps. 
(i) On a representative set of initial data, two differential distributions (his- 
tograms) of the orbits in the computed value of log 10 T L are constructed, 
using two different time intervals for the integration, (m) If the phase space 
of motion is divided [12] , each of these distributions has at least two peaks. 
The peak that shifts (moves in the positive direction in the horizontal axis), 
when the integration time interval is increased, corresponds to the regular 
orbits. The fixed peak (or peaks) corresponds (correspond) to the chaotic 
orbits, (in) The value of log 10 TL at the histogram minimum between the 
peaks is identified, giving a numerical criterion for separating the regular 
and chaotic orbits, (iv) The obtained criterion can be used in any further 
computations to separate the regular and chaotic orbits on much finer initial 
data grids and, rather often, on smaller time intervals of integration. 

We consider two representative sets of initial data. They are formed 
as fine grids in rectangular areas in the (q, e) plane, namely, in the areas 
(0.5 < q < 1.0, < e < 0.9) and (0.6 < q < 0.8, < e < 0.2). The second 
area, which is smaller, is explored with a higher resolution in q and e, i.e., 
the (q, e) grid is finer. 

In Fig. [TJ the distributions (histograms) / of planetary orbits in the com- 
puted value of log 10 T L for the first (large) area are shown for two choices 
of the computation time. The quantity / is defined as the number of orbits 
with log 10 TL in the interval (log 10 TL, log 10 TL + Alog 10 TL), normalized by 
the total number of orbits in the set; Alog 10 TL is set equal to 0.03. As fol- 
lows from the histograms in Fig. [TJ the critical (separating chaos and order) 
value of log 10 T L can be set equal to 2.5. 
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Figure 1: The histograms of circumbinary orbits in the Kepler-IQ system 
in the Lyapunov time logarithm. The histogram of the orbits computed on 
the time interval of 10 3 yr is shown in red, and that computed on the time 
interval of 10 4 yr is shown in blue. 



The resulting stability diagrams are shown in Figs. |2] and |3j Fig. Eh is 
constructed using the Lyapunov exponent criterion, and Fig. [2fr using the 
Lyapunov exponent plus escape-collision criteria. In Fig. |3l the stability 
diagrams are shown in a higher resolution in q and e. 

A direct inspection of Figs. [2fr and [3fr demonstrates that the Lyapunov 
exponent criterion provides a more clear-cut picture of the chaos-order bor- 
ders, in comparison with using the escape-collision criterion at the same time 
interval of integration. What is more, the escape-collision criterion fails to 
identify orbits with high values of the pericentric distance and eccentricity, 
at least at the specified intervals of computation time. The chaos-order bor- 
ders, revealed by the Lyapunov exponent criterion, apparently possess fractal 
structure conditioned by the orbital resonances (discussed below). 

The location of the planet Kepler-I6b is shown in Figs. [2] and [3] by a red 
dot. As it is clear from the diagrams, the planet turns out to be located 
almost at the edge of the chaos domain, in a hazardous vicinity to it - 
just between the instability "teeth". A direct linear extrapolation of these 
"teeth" to the e = axis shows that they correspond to integer resonances 
between the orbital periods of the planet and the central binary. The two 
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Figure 2: The stability diagrams constructed by the Lyapunov exponent (a) and 
by the Lyapunov exponent plus escape-collision (b) criteria. The regular regions 
are shown in white, chaotic in blue, and those for orbits with escapes and collisions 
(at the adopted time interval) in green. The location of the Kepler-16b planet is 
shown by a red dot. 
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Figure 3: Zooms of the stability diagrams in Fig. [21 in a higher resolution. 
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Figure 4: The time behaviour of the resonant argument a of Kepler-lQb. 

teeth surrounding Kepler-l&o correspond to the 5/1 and 6/1 resonances: at 
these resonances, the orbital periods of the planet and the central binary are 
in the ratios 5/1 and 6/1. The smaller teeth centered between the "integer" 
teeth correspond to half-integer resonances. The tooth pointing at Kepler- 
16b corresponds to the 11/2 resonance. 

Why there is no instability in this resonance, at the location of Kepler- 
16b, whereas the neighbouring "integer" teeth extend down to the e = 
axis? 

Chaotic behaviour, which is often present in the dynamics of celestial bod- 
ies, is usually due to interaction of resonances (as in any Hamiltonian system 
[j~2]). How to distinguish between resonant and non-resonant motions? In 
fact, the observable commensurability between the orbital frequencies never 
happens to be completely exact, at least due to observational errors. To 
solve this problem, a "resonant argument" (synonymously, "resonant phase" 
or "critical argument") is introduced. It is a linear combination of some 
angular variables of a system under consideration. In our notations, the res- 
onant argument for an outer resonance (a resonance between a central binary 
and a circumbinary particle) is defined by the formula [U [5] 

a = (k + q)X p - A;A S - lw p , (2) 
where A s and A p are the mean longitudes of a star (in the central binary) and 
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a particle, respectively, and zu p is the longitude of pericenter of the particle; 
k, q and / are integers; q is the so-called resonance order. In outer resonance, 
the ratio of orbital periods of a planet and the central binary is equal to 
(k + q)/k. 

Usually a mean motion resonance splits in a multiplet of subresonances, 
corresponding to a sequence of values of I. This phenomenon is due to pre- 
cession of particle's orbit pericenter, often taking place when the motion is 
perturbed [6j [7J [5] . Here the perturbation is strong, because the mass pa- 
rameter fx = Mij (Mi + M2) (where Mi > M2 are the stellar masses) is large: 
it is equal to 0.227 pp. Hence the precession of the planetary orbit is strong. 

The theory of resonances splitting in multiplets was developed in [5J [7] 
for the high-order inner mean motion resonances in the motion of asteroids. 
According to this theory, inner mean motion resonance k/(k + q) splits in a 
cluster of q + 1 subresonances with I = 0,1, ... ,q. Analogously, the outer 
mean motion resonances also split. 

The coefficients of the subresonant terms in the expansion of the per- 
turbing function are proportional to the eccentricities of the particle and 
perturber in some powers depending on q; in particular, the coefficients of 
the first and last subresonant terms in the multiplet are proportional to the 
eccentricities of the perturber and particle, respectively, in the power equal 
to q [6], [7J . In the pendulum model of subresonance, its width (characteriz- 
ing the "subresonance strength") in the phase space is proportional to the 
square root of this coefficient [T^J [6] . Thus the value of q governs the most 
important properties of the resonant motion. 

Consider two neighbouring outer integer resonances (k + q)/k and (k + 
q + l)/k. The half-integer resonance between them is (2k + 2q + I) /(2k), and 
its order is 2q + 1. Thus, for a high-order half-integer resonance, the power- 
law indices in the subresonant term coefficients are much greater than the 
indices for the neighbouring integer resonances; consequently, the strengths 
of subresonances are much less and their interaction is much weaker. On 
increasing e, they start to overlap much later than in the neighbouring integer 
cases. This explains why the 11/2 resonance at the eccentricity of Kepler- 
16b is stable, whereas the neighbouring integer resonances 5/1 and 6/1 are 
unstable. 

Concluding, Kepler-lQb survives because it is in or close to the half- 
integer 11/2 orbital resonance with the central binary. In the Solar system, 
this phenomenon is analogous to the survival of Pluto and plutinos, which 
are in the 3/2 outer orbital resonance with Neptune. In the latter case the 
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order of the surviving resonance is much smaller, because the mass parameter 
fi in the case of Neptune and Sun is much smaller and therefore the chaos 
border radically shifts to smaller values of semimajor axis (in units of the 
"central binary" size). 

Let us directly inspect how the resonant argument behaves in time. In 
the case of exact resonance, it should librate. The resonant argument a, 
given by Eq. (|2J) , is shown in Fig. H] for the case of main subresonance of 
resonance 11/2 (k = 2, q = 9, I = 0). The time interval spanned by this 
plot is 10 yr, whereas [TJ [2] the binary period is 41.079 d and the planet 
period is 228.78 d. One can see that the argument librates (the amplitude 
of its variation does not exceed 2ir) on the given timescale, but it is evident 
that it slowly circulates on longer timescales. (Note that any choice of I 
different from I = results in a quicker circulation.) It cannot be excluded 
that in reality, described in a more refined (e.g., non- restricted) model of the 
planetary motion, the resonant argument librates on all timescales. In any 
case, the motion is nearly resonant. 

On inspecting the stability diagrams, several immediate questions arise. 
First of all, why Kepler-I6b is so close to the chaos border? For any planet 
in a circumbinary orbit (i.e., at the "outskirts") one expects that at the 
epoch the planet formation the density gradient in the protoplanetary disk 
was negative with a, and thus the density of planetesimals was maximum at 
the chaos border; therefore, the formation of the planet at the border seems 
natural. The migration from the place of birth could be blocked by capture 
in the 11/2 resonance. 

Then, another question arises: why is there only one resonance cell that 
is "occupied"? Why there is not a lot (or several) planets in resonance 
cells like bees in honeycomb? At least for the cells neighbouring to Kepler- 
16b the answer is straightforward. Again, it concerns resonances (and their 
interaction), though quite different from those considered aboveQ These are 
the first-order orbital resonances (k + l)/k of test particles with the planet. 
On increasing k, these resonances start to overlap at some critical k, because 
their widths do not decrease fast enough. On such grounds Wisdom [9] 
inferred that in the case of small eccentricity (e < 0.15) of the particle orbit 
the critical k is given by 

1 In Chirikov's saying, ". . . the physicist first of all looks what are the resonances present 
in a system and how do they interact with each other" [§]. 
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^overlap ~ 0.51/X p , (3) 

where \i v = M p / (M s + M p ) is the mass parameter, and M s and M p are the 
masses of the star and the planet respectively. Using the third Kepler law, 
one finds that k = ^overlap corresponds to the chaotic zone width 

^ ^overlap ~ 1.3/4%, ( 4 ) 

where a p is the semimajor axis of the planet [TOl S]. The particles with a 
within the interval a p ± Aa over iap move chaotically. Due to encounters with 
the planet they escape from this region sooner or later; in such a way a 
particle-free zone around the planet orbit is formed. 

Using Eq. (J3J) and data on a p and masses from [T], and setting M s = 
M 1 + M 2 , one finds for Kepler-lQb: /i p = 3.56 -10" 4 and Aa overlap ~ 0.095 AU. 
Therefore, at least the two neighbouring resonance cells, namely the 9/2 and 
13/2 ones, are purged by the planet residing in the 11/2 resonance cell, 
because they are within the Aa ~ 0.1 AU distance (see Figs. [2] and [3]). Thus 
Kepler-lQb can be called not only a "resonant survivor" but a "resonant 
destroyer" as well. 

Our basic conclusions are as following, (i) The minimum outer border of 
the chaotic domain in the stability diagram for the Kepler-IG system corre- 
sponds to the semimajor 0.60 AU at the zero eccentricity of a test planet, (ii) 
The representative values of the Lyapunov times for the outer orbits in the 
chaotic domain are ~ 3 yr. (ra) The Lyapunov exponent criterion, applied for 
the construction of the stability diagrams, provides much better resolution 
of the chaos-order borders in the stability diagrams in comparison with the 
escape-collision criterion, (iv) The planet Kepler-lQb turns out to be almost 
at the edge of the chaos domain, in a hazardous vicinity to it — just be- 
tween the instability "teeth" in the space of orbital parameters. Kepler-lQb 
survives, because it is in or close to the half-integer 11/2 orbital resonance 
with the central binary. In the Solar system, its survival is analogous to the 
survival of Pluto and plutinos, which are in the 3/2 orbital resonance with 
Neptune. The order of the "occupied" half-integer resonance grows with the 
mass parameter /x of the perturbing binary (in the case of the Solar system, 
the relevant "binary" is Sun and Neptune), because raising \i shifts the sta- 
bility border outwards, (v) The location of a planet in a circumbinary orbit 
close to the stability limit is natural, because one expects a negative with a 
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density gradient in the protoplanetary disk at its outskirts. The neighbour- 
ing (to Kepler-l&o) resonance cells are vacant because they are purged by 
Kepler-16b, due to overlap of first-order resonances with the planet, (vi) No 
radial migration was possible in the Kepler-lQb evolution since its formation, 
because otherwise it would cross the instability "teeth" (see Fig. [3]) and thus 
would be removed. 
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